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SUMMARY 

The thermocaplllary motion of spherical bubbles present In an unbounded 
liquid with a linear temperature distribution Is analyzed, when the Reynolds 
number and the Marangonl number are large. Previous calculations of the ter- 
minal velocity performed for this parametric range did not take Into complete 
consideration the thermal boundary layer present near the surface of the bub- 
ble. In the present study, a scaling analysis Is presented for this problem. 
The thermal boundary layer Is analyzed by an Integral method. The resulting 
terminal velocity Is lower than the one previously calculated, though It Is of 
the same order of magnitude. 


INTRODUCTION 

In the absence of gravity, bubbles present In a host fluid with a tempera- 
ture gradient will move towards the hotter portion of the fluid. This Is due 
to the shear stress at the Interface that Is generated by the temperature 
Induced surface tension gradient. Many studies, both theoretical and experi- 
mental, have been performed on this phenomenon (refs. 1 to 7) where the objec- 
tive Is to determine the terminal velocity of migration and the shape of the 
bubble. 

Two of the most Important parameters that Influence the flow and the heat 
transfer In this problem are the Reynolds number R 0 and the Marangonl number 
Ma. The Marangonl number Is really a Peclet number. Other parameters (Weber 
number. Capillary number, etc.) control the shape of the bubble (ref. 6). Most 
analyses are restricted to small and unit order R 0 and Ma as they use the 
creeping flow equations or perturbations to It. In reference 6, a solution was 
found for any R 0 , so long as Ma Is small. Crespo and Manuel (ref. 7) have 
calculated the terminal velocity for large Ma, where no restriction Is stated 
on R 0 explicitly. However, a large R 0 Is Implied as a thin flow boundary 
layer Is assumed. The authors In reference 7 use a mechanical energy argument, 
first used by Levlch (ref. 8) In his analysis of the rise of bubbles In a 
liquid In a gravitational field, for large Reynold's numbers. In reference 7, 
the temperature field Is not analyzed at all and the work of the surface ten- 
sion forces at the Interface that was calculated Is suspect, especially at the 
front and rear stagnation points. In the present study, the analysis In 
reference 7 Is extended and the temperature field Is analyzed by an Integral 
method. This results In a decrease In the calculated terminal velocity com- 
pared to the value reported In reference 7. 

NOMENCLATURE 

A temperature gradient far away from the bubble 
Ca Capillary number, viVr/<* 



f reference quantity for the ratio V'/V^ 

p 

Ha Harangonl number, (-Oj) AR.|/(ya) 

Pr Prandtl number, v/a 
p,P dimensionless and dimensional pressure 
bubble radius 

Re Reynolds number, VpR^/v 

2 

R surface tension Reynolds number, (-<f T )AR,/(tiv) 

<3 I I 

r,R dimensionless and dimensional radial coordinate 
T temperature 

T‘ dimensionless transformed temperature 
T s bubble surface temperature 

t time 

u, U dimensionless and dimensional radial velocity 

v, V dimensionless and dimensional tangential velocity 

V velocity vector 

2 

Wb Weber number, />V R R^ /«? 

a thermal dlffuslvlty 

& reference quantity for the boundary layer thickness 

e tangential coordinate 

x Inverse of the boundary layer thickness at any e 
v viscosity 

v kinematic viscosity 

S stretched radial coordinate 

p density 

a surface tension 

aj temperature coefficient of surface tension 

t viscous stress 

<t> azimuthal coordinate 
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Superscripts 


' correction fields In the boundary layer 

Subscripts 

1 Invlscld fields 

R reference values 

» values far from the bubble 


Formulation 

As the bubble moves through the fluid at Its terminal velocity. It Is con- 
venient to choose a coordinate system on the bubble with the origin at Its 
center of mass. In this system, the bubble Is stationary and the fluid outside 
approaches the bubble at the terminal velocity. The velocity field Is steady. 
The temperature field Is not steady, as the bubble constantly moves to a warmer 
region. However, the gradient of the temperature field Is steady and It Is 
easy to transform to another temperature that Is steady (ref. 5). 


The flow Is considered to be Incompressible and laminar. All physical 
properties other than the surface tension are taken to be Independent of tem- 
perature and are hence spatially constant. The bubble Is assumed to retain Its 
spherical shape. From the geometry and the boundary conditions, the problem 
Is symmetric about the flow direction. It Is assumed that the viscosity and 
thermal conductivity of the gas Inside the bubble are negligible compared to 
those In the liquid outside It. Hence, the flow and the heat transfer within 
the bubble Is not analyzed. The coordinate system Is R,e,$ with the origin 
at the center of the bubble (fig. 1). © Is measured counterclockwise from the 

point of Incidence of the flow. The basic equations and boundary conditions 
describing the flow are the same as In reference 6 and are reproduced below In 
equations (1) to (10). 
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where R-j Is the radius of the bubble, A Is the temperature gradient In the 
liquid far away from the bubble, Vr = (-<j-r)AR-|/y, Is a reference velocity 
determined from the shear stress condition at the Interface, T* Is the 
dimensionless transformed temperature which Is steady. _In what follows, 
unless otherwise mentioned, T will be used to denote T‘, for convenience. 

<*T = do/dT Is the rate of change In surface tension with temperature and Is 
taken to be a constant. It Is usually negative. Vo, Is the terminal velo- 
city. The basic equations are 


1 

r 


_ i_ 
2 ar 


(r 2 u) 


+ TliTe ae < v s1n 6 > 


3 


0 


( 2 ) 


(3) 


2 

au . y au _ v_ 
dr r 30 r 


u 



v 3v uv 

r ae r 


+ I_ [ V 2 U _ 2 

ar R o L 


2y 

2 


2_ av 

_2 ae 


lie , 
"r ae 


<* L 


V 2 V + — — 

r 2 ae 


■ ? •] 


r 2 s 


v + u 



v aT 
r ae 


l_ 

Ha 


T 




1 

r 2 s1n e 


a_ 

ae 



(4) 


(5) 


where R 0 * VpR-|/v * (-<jj)AR?/(vv) Is the surface-tension Reynolds number, 
Ma * Pr R 0 Is the Harangonl number, and v® * V®/Vr. 

The boundary conditions are 


At r = 1, 


As r -» ®, 

u -» -v cos e, 
00 ’ 


u ■ 0 

a / v\ aT 
ar \r) 55 ae 



v -» v sin e, 

CO 



T •* r cos e 


( 6 ) 

(7) 

( 8 ) 


(9) 


In addition to the above conditions, we have another condition resulting 
from the fact that the total force acting on the bubble Is zero, as the bubble 
Is moving at a constant velocity. This condition Is used to determine v®, 
the dimensionless terminal velocity. 


Jf [’Re 5 '" 2 o - t rr cos 9 s,n °] R _ Ri 


de = 0 


( 10 ) 


where (x Re ) 


R-R, ‘ “[" » B (") 


and < t rr) 

R=R 1 KK R=R 1 


= 2y 


/au\ 

V 3 R /r=r, 


P(R r e) 


As an alternative to the net force condition, Crespo and Manuel (ref. 7), 
following Levlch (ref. 8), have used a condition arising from the conservation 
of mechanical energy. In a reference frame In which the bubble Is stationary, 
conservation of mechanical energy Implies that the rate of energy dissipated 
In the liquid by viscous forces must equal the rate of work done by the surface 
tension forces at the Interface. Hence, this condition may be written as 
(refs. 7 to 9). 
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Simplifying and using U = 

£ ( v SL, s,n 9 de 



(V x curl V) - \ J- (V 
this becomes 

(RV) - |^| 2 sin © dR de 



ds 


+ 2 V R^ (V 2 ) R=R ^ sin e de 


( 12 ) 


2 2 

where [a/aR(RV) - au/ae] = | cur 1 V| = the square of the magnitude of the 
vortlclty and T Is the dlmenslonaT transformed temperature. 


While In typical boundary layer problems (like flow over a flat plate), 
using an energy method such as the one above to determine the drag Is more 
difficult than directly Integrating the surface stresses. It appears that In 
problems with spherical bubbles, the energy method Is easier to apply. Levlch 
(ref. 8) and Moore (ref. 10) have used such a method to predict the terminal 
velocity of rise of bubbles In a gravitational field, when the Reynolds number 
Is large. 


Analysis 

The boundary value problem that has been formulated Is very difficult to 
solve. In the Introduction, the various conditions have been mentioned under 
which solutions have been obtained. In the present study, our Interest Is to 
predict the terminal velocity for large R a and Ma. As a first step, an 
estimate for the terminal velocity, l.e., Vr Is determined for large R a and 
Ma. The previous estimate Vr * ( -oj) AR^ /p (ref. 6) Is valid only for small 
Ma, as the shear stress was assumed to be of order uVr/R-| , which excludes the 
presence of sharp gradients, l.e., the presence of boundary layers. We will 
determine Vr by scaling analysis. For large R<, and Ma, we expect a 
boundary layer to be present near the bubble surface. We will represent the 
velocity field V as the sum of the Invlscld velocity field V^ and a boun- 
dary layer correction velocity field V', l.e., 

i-i H3) 

The order of magnitude of V^ Is denoted to be V R and that of V' to 
be f V r . The correction field V' Is nonzero within the boundary layer and Is 
zero outside It. 6 denotes thickness of the boundary layer. Vr, f, and 6 
are all unknown and must be determined. These are determined below by appro- 
priate balance of terms In the equations and boundary conditions. 
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Nondlmenslonal quantities are defined as follows 
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where U^, V^, and are the dimensional Invlscld velocity components and 
the pressure and U 1 and V' are the correction velocities In the boundary 
layer. Using equation (14) In equations (2) to (5), the following equations 
are obtained 
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The shear stress condition (eq. (7)) becomes 

_ ( - tf T )AR 1 31 
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f 3 _ / V 1 \ 

* 35 V + * 5 A =0 " 
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The mechanical energy condition (eq. (12)) becomes 
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( 20 ) 
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Three conditions are needed to determine the unknowns 6, f, and Vr. 
These are described below 

(1) Inertia and viscous forces are of the same order of magnitude In the 

, ,, . u 3v 1 1 1 

boundary layer. Hence u^ — 


(1 ♦ *5) 


2 35 


’ __.2 3 v 1 ~ 1 

(1 + 65) ^ J. 


Since u 1 ~ 6 for 5 ~ 1 , * ~ * « ~ 

(2) The left-hand side of the mechanical energy condition (eq. (21)) 
represents the rate of work due to surface tension forces. The first term on 
the right-hand side represents the energy dissipated In the boundary layer 
(since curl V * 0 outside the boundary layer) and the second term represents 
the energy dissipated by the potential flow and an additional boundary layer 
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contribution (note that the flow considered here Is potential because It Is 
Irrotatlonal and not because It Is Invlscld; hence the potential flow does 
dissipate energy). We expect f ~ 0(1) or less and i « 1. Hence f 2 /« 
represents the largest portions of the energy dissipated In the boundary layer. 
Hence 

(-*T )AR 1 ~ larger of I and 1 . 

wV r L* -I 

(3) Similarly, from the shear stress condition, 

f r ( -"T )AR l 1 

larger of I — ^ and 1J 

It Is assumed that (dT/de) r _i ~ 1. Physically, this means that the tem- 
perature between the two stagnation points on the bubble surface Is of order 
AR-j . Solving these three conditions for 6, f, and Vr, we find that there 
are two possibilities 

a) J— 1 < 22 > 


Since & must be small compared to one. 


<I» ‘ - ^73 • 


Again, 6 « 1 => R„ » 1 . 

For the first possibility, f 2 /« * l/V^o which Is small compared to 
one. Hence all the work that Is done by the surface tension forces Is dissi- 
pated by the potential flow, with negligible amounts dissipated In the boundary 
layer. For the second possibility, the converse Is true and hence the predom- 
inant dissipation of energy occurs within the boundary layer, with negligible 
amounts dissipated outside It by the potential flow. Also, according to the 
first possibility, the correction velocity within the boundary layer Is small 
compared to the potential flow velocity, whereas, according to the second pos- 
sibility, the two are comparable. Crespo and Manuel assumed that possibility I 
Is correct, as they neglected energy dissipation In the boundary layer, refer- 
ring to arguments by Levlch (ref. 8, Ch 7, § 82) for justification. The velo- 
city scale In the second possibility Is the one frequently used In literature 
for thermocaplllary flows with large R„ and was first obtained by Ostrach 
(ref. 11). 


R » 1. 
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V R 


(-Oj^A^v 
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(23) 


So far, the author has not been able to produce an argument which tells a 
priori which of the two possibilities Is the one that Is physically correct and 
Is chosen by nature. Therefore, both possibilities were pursued. It was found 
that the second possibility, analyzed by a boundary layer Integral method simi- 
lar to the Karman-Pohlhausen method, was never able to predict any terminal 
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velocity. While the scaling analysis for possibility II revealed that all the 
work done by the surface tension forces could be dissipated within the boundary 
layer, calculations by the Integral method revealed that for any terminal velo- 
city, the boundary layer dissipation was always less than (and hence could 
never equal) the rate of work done by the surface tension forces. Thus possi- 
bility II was rejected. Therefore, In what follows, possibility I Is chosen 
to be physically correct and will be pursued In the remainder of this work. 


Using equation (22) In the energy equation (eq. (18)) and the mechanical 
energy condition (eq. (21)), we may neglect terms containing u' and v', for 
large R tf . The following simplified problem for T and Is then 

obtained 


v + u 31 + M ai _ 1. 2 
» U 1 dr + r ae Ma v 
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As r -» ob. 
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ar 
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sin e de 



< v 5> 


r*l 


sin e de 


-v cos 
00 




a v sin e 
00 


(- 



(27) 


(28) 


Solution 

The above equation for T, equation (24), though linear. Is not easy to 
solve. Crespo and Manuel (ref. 7) neglected the right-hand side of 
equation (24) for large Ma, yielding 


v 

CO 


+ u, 


31 

3r 


_i II 
r ae 


0 


Evaluating this equation at r * 1, they obtained 



se 


-V 


oo 


Using this In equation (27), the terminal velocity v Is obtained as 

00 


(29) 


(30) 


v 


1 

3 


(31) 


Thus, they were able to calculate the terminal velocity without solving for the 
temperature field. However, equation (29) Is not strictly valid as It cannot 
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satisfy equation (25) at the bubble surface. Therefore, a thermal boundary 
layer must be present near the bubble surface. Also equation (30) Is Incor- 
rect, especially at the stagnation points, where the left-hand side Is zero. 
We will solve equation (24) by an approximate Integral method. Equation (24) 
may be Integrated with respect to r to yield 



[r(v^ + «v')T sin 6 


2 2 3 

v^r cos e sin e] dr + ^- (cos e - 


1 _ 

Ha 


(sine 


II 

ae 


+ r 


s1n 2 e^dr | 


(32) 


Evaluating this equation at e « 0 and assuming that the Integrals are zero, 
Ci * 1. However, at e = *, the second term Is nonvanishing. Hence, one 
of the Integrals must be finite at e = *. Physically, there Is a thermal 
wake behind the bubble, for, as the bubble moves Into warmer regions. It dis- 
places the warmer fluid. The displaced energy has to be convected downstream 
In a thermal wake. The Integral In the left-hand side of equation (32) must 
be nonzero at e « *, to account for the energy In the wake. 


The following temperature distribution Is assumed 
T(r,6) - r cos e - j + ^3 (X - 3) [ XT s (0) " * 2 ^ cos e ] 

♦ (x z 3) |f cos 6 - 3T s (e)je _X(0)(r " 1) (33) 

where T s (e) « T(l,e) Is the transformed steady temperature at the bubble 
surface. T s and x are unknowns to be determined. The assumed temperature 
distribution satisfies the boundary conditions (eqs. (25) and (26)) and makes 
the Integrals In equation (32) finite. The exponential function represents a 
thermal boundary layer (X Is the Inverse of the boundary layer thickness) and 
also accommodates the energy In the wake. The remaining terms In equation (33) 
represent the temperature distribution outside the boundary layer. 


Equations (32) and (24) evaluated at r * 1 are used to determine T s 
and X. 


8(X - 3) 


U s - cos * 11 [\ 

J X U - 3) ' 


cos 0 - 3T 


) 


(j »s » - 3T S ) 1 1 | (C °|^ - J) ■ 0 <34 > 
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( 35 ) 


v . ♦ 2 *. s1 " 6 T i 


k{- 


3 cos e + 


12 


(X - 3) 


t 


XT. - 


(X + 4) 


cos e 


■] 


[! 


cos e - 3T S * T ; * cot 9 T $ 


.] 


T'(0) = T'(ir) = 0 from symmetry conditions 


‘! 


(36) 


where T$ = dT s /de and T$ = d 2 T s /de 2 . 1/Ma terms In equation (34) have been 
neglected. They are retained In equation (35), as without them, the equation 
Is singular at 6=0 and w. This system Is still not easy to solve. Hence, 
we will get a solution by satisfying the equations only at 6=0 and ir/2. 


It Is difficult to satisfy the equations at 
and will not be attempted. 


6 = *, as T s Is steep there 


Is chosen to be 


T s (6) 


a + b 


(=) ! • • (;)' 


(37) 


This satisfies the condition T s (0) = 0. The six unknowns a, b, c, x(0), 
x(*/2) and v„ are solved from equations (34) and (35) evaluated at 6=0 


and ir/2, the condition 
(eq. (27)). 


Tc(*) = 0 and the mechanical energy condition 


The solution for large Ha Is 


T s (0) - a - ff 


. 3c 13 

b - - r = - r 


X(0) = 0.53 Ma 


(?)■•■ 


009 Ma 


v = 


_ 1 3(» + 3) 

12m 2 


0.235 


AT $ = T s (0) - T $ (ir) = 1.44 


1 


(38) 


RESULTS AN0 DISCUSSION 

The results for the terminal velocity, viz, = 0.235 (-«jj)ARi/u Is 
lower than that obtained by Crespo and Manuel, viz, V* = l/3((-<rf)ARi/v) . How- 
ever, even though their approach Is not fully justified, both results are of 
the same order of magnitude and only differ In the multiplicative constant In 
the expression for The main difference Is that In the present study, the 

thermal boundary layer has been treated more completely In arriving at the 
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result for the terminal velocity. For large Ma, v m Is Independent of the 
Marangonl number. The result for v a also supports the scaling analysis that 
has been performed, as the estimate for v* Is modified by a multiplicative 
constant that Is only of unit order of magnitude. The temperature on the 
bubble surface at the front stagnation point Is T s (0) = 1.01 and Is only 
slightly different from the temperature of the free stream (T = 1) at the same 
axial location. The thickness of the thermal boundary layer varies as 1/Ha. 
The boundary layer has a small thickness (1.89/Ma) at 0=0, Is about 58 times 
thicker at e = */2 and Is Infinitely thick (from eq. (34)) at 0 = *. An 
Interesting conclusion that Is supported by the present analysis Is the fact 
that Ma Is a singular perturbation parameter for this problem, which has been 
recognized and considered by Subramanlan (ref. 5). For Ha = 0, the solution 
to T Is (ref. 6) 


T(r,0) = r cos 0 + ^ (39) 


Comparing this to T(r,0) In equation (33), we see that the sign of the second 
term Is reversed. The coefficient -1/2 for this term In equation (33) was 
determined by requiring that the Integral In the left-hand side of equation (32) 
be finite, l.e., that the motion of the bubble does not create an Infinite flux 
of convected energy. Hence, the coefficient of this term must be -1/2 for 
any nonzero Ma, as otherwise, the flux of energy convected would be Infinite. 
Since the coefficient Is +1/2 for Ma = 0, we conclude that this problem Is 
singular with respect to perturbations In Ma (l.e., the Inclusion of energy 
convection terms), as the presence of convection drastically changes the nature 
of the temperature distribution as r 

For small Marangonl, Weber, and Capillary numbers and small W«, the 
shape of the bubble was obtained In reference 6 to be 

1 c vVp ? 

n (e) - - (3 cos^0 - 1) (40) 


where the bubble surface Is located at r = 1 + n. The shape of the bubble 
represented by equation (40) Is a spheroid with Its minor axis In the flow 
direction. The same result Is also expected to be valid for the large R d 
and Ma that Is being considered In the present study for small Wb, Ca, and 
Wo, because the flow boundary layer Is thin and Introduces only small 
changes to the velocity and pressure fields, compared to these fields In poten- 
tial flow. Since potential flow fields were used In reference 6 to obtain 
equation (40), the shape of the bubble for the two cases must be the same. 

This result Is not expected to be valid In the vicinity of the rear stagnation 
point of the bubble, as we expect the boundary layer thickness to be Infinite 
there. 
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Figure 1. - Sketch of a migrating droplet. 
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